Our goal is to approximate a bitmap using a smooth function \( u \)
\( \color{red}{u(x,y)} \approx BITMAP(x,y) \).
Our slim feature \( u \) it is formed by the linear combination of many base functions
\( u(x,y) = \sum_{i=1,j=1}^{N_x,N_y} u_{i,j} B^x_{i}(x) B^y_{j}(y) \) where \( u(x,y) = _{j}(y) \) are the coefficients of this combination, the values which we must calculate.
Our base functions are spread over \( N_x*N_y \) mesh elements. So there are also \( N_x*N_y \) base functions, and therefore we need to calculate \( N_x*N_y \) coefficients \( u_{i,j} \). Each of these coefficients \( u_{i,j} \) tells me how to scale one base function, a two dimensional B-spline function \( B_{i,j;2}(x,y)=B^x_{i}(x) B^y_{j}(y) \).
How do we calculate these \( \{u_{i,j}\}i=1,..,N_x;j=1,...,N_y \) coefficients?.
Instead of choosing pixel values, we will use the averaging method. We also use our B-spline base functions for averaging. We want our approximation to approximate the bitmap \( u(x,y) \approx BITMAP(x,y) \).
So let's choose any B-spline function
\( v(x,y)=B_{k,l;2}(x,y)=B^x_{k}B^y_{l} \). We chose a function with indices \( k,l \). For simplicity, we name it \( v \).
Three exemplary selections of our B-spline function are illustrated in the figure. We will call this function a test function. This is the standard nomenclature used in the finite element method.
We then use our test function to compute the average of the bitmap surrounding the function.
In other words, we want our approximation \( u \) to approximate our bitmap with such averaging distribution as our test function gives us. To write this operation mathematically, we multiply our approximation equations by the test function and calculate the integral
\( \int_{\Omega} u(x,y) v(x,y) = \int_{\Omega} BITMAP(x,y) v(x,y) \).
The focal point has the greatest impact on the value of this mean, according to the appearance of the B-spline basis function, and the more distant points have a smaller effect on the mean value.
Now we can note that one selection of the test function
\( v \) gives us one equation.
How many ways can we choose our test functions? Since each of them is a B-spline associated with the center of an element, and there are \( N_x*N_y \) elements, it means we can choose \( N_x*N_y \) test functions.
Each of these selections of test functions gives us one equation, so we get a system of \( N_x*N_y \) equations
\( \int u(x,y) \color{blue}{B^x_1(x)B^y_1(y)} = \int BITMAP(x,y) \color{blue}{B^x_1(x)B^y_1(y)} \)
\( \int u(x,y) \color{blue}{B^x_1(x)B^y_2(y) } = \int BITMAP(x,y) \color{blue}{B^x_1(x)B^y_2(y)} \)
\( \vdots \)
\( \int u(x,y) \color{blue}{B^x_k(x)B^y_l(y)} = \int BITMAP(x,y) \color{blue}{B^x_k(x)B^y_l(y)} \)
\( \vdots \)
\( \int u(x,y) \color{blue}{B^x_{N_x}(x)B^y_{N_y-1}(y)} = \int BITMAP(x,y) \color{blue}{B^x_{N_x}(x)B^y_{N_y-1}(y)} \)
\( \int u(x,y) \color{blue}{B^x_{N_x}(x)B^y_{N_y}(y)} = \int BITMAP(x,y) \color{blue}{B^x_{N_x}(x)B^y_{N_y}(y)}. \)
The question that needs to be asked now is that of the unknowns in our system of equations. We have \( N_x*N_y \) equations, but where are the unknowns? To reach the unknowns, remember how our function is constructed \( u \) approximating a bitmap.
\( \color{red}{u(x,y) \approx \sum_{i,j} u_{i,j} B^x_{i}(x) B^y_{j}(y)} \)
Our unknowns are coefficients \( \{u_{i,j}\}i=1,...,N_x,j=1,...,N_y \)
specifying how we "pull" individual B-splines "hills" in order to obtain the best bitmap approximation.
So we put our linear combination of B-spline in place \( u \) and we get

\( \int \color{red}{\sum_{i,j} u_{i,j} B^x_{i}(x) B^y_{j}(y)} \color{blue}{B^x_1(x)*B^y_1(y)} = \int BITMAP(x,y) \color{blue}{B^x_1(x)*B^y_1(y)} \)

\( \int \color{red}{\sum_{i,j} u_{i,j} B^x_{i}(x) B^y_{j}(y)} \color{blue}{B^x_1(x)*B^y_2(y)} = \int BITMAP(x,y) \color{blue}{B^x_1(x)*B^y_2(y)} \)

\( \vdots \)

\( \int \color{red}{\sum_{i,j} u_{i,j} B^x_{i}(x) B^y_{j}(y) } \color{blue}{B^x_k(x)*B^y_l(y)} = \int BITMAP(x,y) \color{blue}{B^x_k(x)*B^y_l(y) } \)

\( \vdots \)

\( \int \color{red}{\sum_{i,j} u_{i,j} B^x_{i}(x) B^y_{j}(y)} \color{blue}{B^x_{N_x}(x)*B^y_{N_y-1}(y)} = \int BITMAP(x,y) \color{blue}{B^x_{N_x}(x)*B^y_{N_y-1}(y)} \)

\( \int \color{red}{\sum_{i,j} u_{i,j} B^x_{i}(x) B^y_{j}(y)} \color{blue}{B^x_{N_x}(x)*B^y_{N_y}(y)} =\int BITMAP(x,y) \color{blue}{B^x_{N_x}(x)*B^y_{N_y}(y)} \)

We now use the algebraic property of the system of equations, which allows us to extract the sum with the coefficients over the individual integrals
\( \color{red}{ \sum_{i,j} u_{i,j}} \int \color{red}{B^x_{i}(x) B^y_{j}(y)} \color{blue}{B^x_1(x)*B^y_1(y)} = \int BITMAP(x,y) \color{blue}{B^x_1(x)*B^y_1(y) } \)

\( \color{red}{\sum_{i,j} u_{i,j}} \int \color{red}{B^x_{i}(x) B^y_{j}(y)} \color{blue}{B^x_1(x)*B^y_2(y) } = \int BITMAP(x,y) \color{blue}{B^x_1(x)*B^y_2(y) } \)

\( \vdots \)

\( \color{red}{\sum_{i,j} u_{i,j}} \int \color{red}{B^x_{i}(x) B^y_{j}(y)} \color{blue}{B^x_k(x)*B^y_l(y) } = \int BITMAP(x,y) \color{blue}{B^x_k(x)*B^y_l(y)} \)

\( \vdots \)

\( \color{red}{\sum_{i,j} u_{i,j}} \int \color{red}{B^x_{i}(x) B^y_{j}(y)} \color{blue}{B^x_{N_x}(x)*B^y_{N_y-1}(y)} = \int BITMAP(x,y) \color{blue}{B^x_{N_x}(x)*B^y_{N_y-1}(y)} \)

\( \color{red}{\sum_{i,j} u_{i,j}} \int \color{red}{B^x_{i}(x) B^y_{j}(y)} \color{blue}{B^x_{N_x}(x)*B^y_{N_y}(y)} =\int BITMAP(x,y) \color{blue}{B^x_{N_x}(x)*B^y_{N_y}(y)} \)

Note now that the above system of equations can be written in matrix form.
In particular, coefficients marked in red \( \color{red}{u_i,j} \) form a vector of unknowns, \( \int \color{red}{B^x_{i}(x) B^y_{j}(y)}\color{blue}{B^x_k(x)*B^y_l(y)} dx \) integrals form a matrix whose lines are numbered \( i=1,...,N_x; j=1,...,N_y \) , and columns from \( k=1,...,N_x; l=1,...,N_y \). This is because the integrals contain the products of two two-dimensional B-spline functions, two "hills" spanning a finite element mesh. The first B-spline function, highlighted in black, is used for approximation, while the second B-spline function, highlighted in blue, is used to sample the bitmap (to compute the average of the bitmap). The vector of the right-hand side is the vector of integrals in which the bitmap is multiplied by the "blue" testing B-spline used for averaging, and decomposed according to the distribution described by our test B-spline.
In summary, we have thus obtained the following system of equations in which the unknowns are the coefficients
\( \color{red}{\{u_{i,j}\}i=1,...,N_x;j=1,...N_y} \)
\( \begin{bmatrix} \int{\color{red}{B^x_{1,p}B^y_{1,p}}}\color{blue}{B^x_{1,p}B^y_{1,p}} & \int{\color{red}{B^x_{1,p}B^y_{1,p}}}\color{blue}{B^x_{2,p}B^y_{1,p}} & \cdots & \int{\color{red}{B^x_{1,p}B^y_{1,p}}}\color{blue}{B^x_{N_x,p}B^y_{N_y,p}} \\ \int{\color{red}{B^x_{2,p}B^y_{1,p}}}\color{blue}{B^x_{1,p}B^y_{1,p}} & \int{\color{red}{B^x_{2,p}B^y_{1,p}}}\color{blue}{B^x_{2,p}B^y_{1,p}} & \cdots & \int{\color{red}{B^x_{2,p}B^y_{1,p}}}\color{blue}{B^x_{N_x,p}B^y_{N_y,p}} \\ \vdots & \vdots & \vdots & \vdots \\ \int{\color{red}{B^x_{N_x,p}B^y_{N_y,p}}}\color{blue}{B^x_{1,p}B^y_{1,p}} & \int{\color{red}{B^x_{N_x,p}B^y_{N_y,p}}}\color{blue}{B^x_{2,p}B^y_{1,p}} & \cdots & \int{\color{red}{B^x_{N_x,p}B^y_{N_y,p}}}\color{blue}{B^x_{N_x,p}B^y_{N_y,p}} \\ \end{bmatrix} \begin{bmatrix} \color{red}{u_{1,1}} \\ \color{red}{u_{2,1}} \\ \color{red}{\vdots} \\ \color{red}{u_{N_x,N_y}}\\ \end{bmatrix} \)
\( = \begin{bmatrix} \int BITMAP(x,y) \color{blue}{B^x_1(x)*B^y_1(y)} dx \\ \int BITMAP(x,y) \color{blue}{B^x_1(x)*B^y_2(y)} dx \\ \vdots \\ \int BITMAP(x,y) \color{blue}{B^x_{N_x}(x)*B^y_{N_y}(y) } dx \\ \end{bmatrix} \)
In a "magical" way, solving the above system of equations will give us the best approximation factors for the bitmap through our B-splines.
The aforementioned "magic" results from determining the individual equations of our system, in which we postulated that on average - in the sense of the distribution determined by the B-spline testing functions - our approximation will correspond to the bitmap.